# 确定维数

# step9 Determine the 'dimensionality' of the dataset
# NOTE: This process can take a long time for big datasets, comment out for expediency. More
# approximate techniques such as those implemented in ElbowPlot() can be used to reduce
# computation time

pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)


ElbowPlot(pbmc)


jackstraw [ˈdʒækstrɔː] n. 稻草人；小木片（游戏用）




() ElbowPlot()
./seurat-4.1.0/R/visualization.R:3667:ElbowPlot <- function(object, ndims = 20, reduction = 'pca') {

选择相关的维度，这个比 Jackstraw 要快得多。

#' Quickly Pick Relevant Dimensions
#'
#' Plots the standard deviations (or approximate singular values if running PCAFast)
#' of the principle components for easy identification of an elbow in the graph.
#' This elbow often corresponds well with the significant dims and is much faster to run than
#' Jackstraw
#'
#' @param object Seurat object
#' @param ndims Number of dimensions to plot standard deviation for
#' @param reduction Reduction technique to plot standard deviation for
#'
#' @return A ggplot object
#'
#' @importFrom cowplot theme_cowplot
#' @importFrom ggplot2 ggplot aes_string geom_point labs element_line
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' ElbowPlot(object = pbmc_small)
#'
ElbowPlot <- function(object, ndims = 20, reduction = 'pca') {
  # 获取每个PC对应的标准差
  data.use <- Stdev(object = object, reduction = reduction)
  # 如果长度为0，报错
  if (length(x = data.use) == 0) {
    stop(paste("No standard deviation info stored for", reduction))
  }

  # 如果需要显示的维度 大于 标准差的个数，则警告，让ndim为现有维度
  if (ndims > length(x = data.use)) {
    warning("The object only has information for ", length(x = data.use), " reductions")
    ndims <- length(x = data.use)
  }

  # y label
  stdev <- 'Standard Deviation'

  # 画图，拼凑数据框, dims, stdev 两列
  plot <- ggplot(data = data.frame(dims = 1:ndims, stdev = data.use[1:ndims])) +
    # x是1,2,... ， y是每个PC对应的方差
    geom_point(mapping = aes_string(x = 'dims', y = 'stdev')) +
    # x轴标签 是去掉下换线后缀后的key
    labs(
      x = gsub(
        pattern = '_$',
        replacement = '',
        x = Key(object = object[[reduction]])
      ),
      y = stdev
    ) +
    theme_cowplot()
  # 返回结果
  return(plot)
}








()Stdev()
./seurat-object-4.0.4/R/generics.R:755:Stdev <- function(object, ...) {

#' Get the standard deviations for an object
#'
#' @param object An object
#' @param ... Arguments passed to other methods
#'
#' @return The standard deviations
#'
#' @rdname Stdev
#' @export Stdev
#'
#' @concept data-access
#'
Stdev <- function(object, ...) {
  UseMethod(generic = 'Stdev', object = object)
}
普通的S3泛型函数。



./seurat-object-4.0.4/R/seurat.R:1726:Stdev.Seurat <- function(object, reduction = 'pca', ...) {

#' @param reduction Name of reduction to use
#'
#' @rdname Stdev
#' @export
#' @method Stdev Seurat
#'
#' @examples
#' # Get the standard deviations for each PC from the Seurat object
#' Stdev(object = pbmc_small, reduction = "pca")
#'
Stdev.Seurat <- function(object, reduction = 'pca', ...) {
  CheckDots(...)
  return(Stdev(object = object[[reduction]]))
}
获取内部的对象，然后抛给他们干活。

> pbmc_small[['pca']]
A dimensional reduction object with key PC_ 
 Number of dimensions: 19 
 Projected dimensional reduction calculated:  TRUE 
 Jackstraw run: TRUE 
 Computed using assay: RNA



./seurat-object-4.0.4/R/dimreduc.R:402:Stdev.DimReduc <- function(object, ...) {

#' @rdname Stdev
#' @export
#' @method Stdev DimReduc
#'
#' @examples
#' # Get the standard deviations for each PC from the DimReduc object
#' Stdev(object = pbmc_small[["pca"]])
#'
Stdev.DimReduc <- function(object, ...) {
  CheckDots(...)
  return(slot(object = object, name = 'stdev'))
}

该函数包装好几层，最后就是返回个 dimReduc@stdev

> (pbmc_small[['pca']])@stdev
 [1] 2.7868782 1.6145733 1.3162945 1.1241143 1.0347596 0.9876531 0.8501773 0.8225654 0.7607182 0.7387425
[11] 0.6680697 0.6216919 0.5463279 0.4933780 0.4639999 0.4404766 0.3246018 0.2762482 0.1704923














()JackStraw
#pbmc <- JackStraw(pbmc, num.replicate = 100)
./seurat-4.1.0/R/dimensional_reduction.R:49:JackStraw <- function(

确定每个PC的显著程度。
- 随机抽取 HVG 基因子集
- 计算 这些随机基因的 projected PCA scores
- 比较 随机基因的 PCA scores 和 观测到的PCA scores，决定显著性。
- 结果是 每个基因和每个PC的相关性的p值


#' Determine statistical significance of PCA scores.
#'
#' Randomly permutes a subset of data, and calculates projected PCA scores for
#' these 'random' genes. Then compares the PCA scores for the 'random' genes
#' with the observed PCA scores to determine statistical signifance. End result
#' is a p-value for each gene's association with each principal component.
#'
#' @param object Seurat object
#' @param reduction DimReduc to use. ONLY PCA CURRENTLY SUPPORTED.
#' @param assay Assay used to calculate reduction.
#' @param dims Number of PCs to compute significance for
#' @param num.replicate Number of replicate samplings to perform #抽样次数
#' @param prop.freq Proportion of the data to randomly permute for each replicate #每次抽样百分比
#' @param verbose Print progress bar showing the number of replicates # 是否打印重复次数的进度
#' that have been processed.
#' @param maxit maximum number of iterations to be performed by the irlba function of RunPCA #计算PCA的irlba()最大迭代次数
#'
#' @return Returns a Seurat object where JS(object = object[['pca']], slot = 'empirical')
#' represents p-values for each gene in the PCA analysis. If ProjectPCA is
#' subsequently run, JS(object = object[['pca']], slot = 'full') then
#' represents p-values for all genes.
#'
#' @importFrom methods new
#' @importFrom pbapply pblapply pbsapply
#' @importFrom future.apply future_lapply future_sapply
#' @importFrom future nbrOfWorkers
#'
#' @references Inspired by Chung et al, Bioinformatics (2014)
#' @concept dimensional_reduction
#'
#' @export
#'
#' @examples
#' \dontrun{
#' data("pbmc_small")
#' pbmc_small = suppressWarnings(JackStraw(pbmc_small))
#' head(JS(object = pbmc_small[['pca']], slot = 'empirical'))
#' }
#'
JackStraw <- function(
  object,
  reduction = "pca",
  assay = NULL,
  dims = 20,
  num.replicate = 100, #默认重复100次
  prop.freq = 0.01, #默认抽样 1% 的细胞
  verbose = TRUE,
  maxit = 1000
) {
  # (A1) 现在只支持PCA降维
  if (reduction != "pca") {
    stop("Only pca for reduction is currently supported")
  }

  # (A2) 根据参数指定 apply 函数版本
  # verbos 且 默认单核时，显示进度条
  if (verbose && nbrOfWorkers() == 1) {
  	# pbapply: Adding Progress Bar to '*apply' Functions
    my.lapply <- pblapply
    my.sapply <- pbsapply
  } else {
  	# 否则，使用多线程版本的 apply 函数，见 补充说明1
    my.lapply <- future_lapply
    my.sapply <- future_sapply
  }

  # (A3) assay 如果为空，则使用默认assay
  assay <- assay %||% DefaultAssay(object = object)

  #(A4) 如果是 sctransform 处理的，则报错。
  if (IsSCT(assay = object[[assay]])) {
    stop("JackStraw cannot be run on SCTransform-normalized data.
         Please supply a non-SCT assay.")
  }

  # (A5) 显示的维度，不能超过实际PC个数
  if (dims > length(x = object[[reduction]])) {
    dims <- length(x = object[[reduction]])
    warning("Number of dimensions specified is greater than those available. Setting dims to ", dims, " and continuing", immediate. = TRUE)
  }

  #(A6) 显示的维度，不能超过细胞数
  if (dims > nrow(x = object)) {
    dims <- nrow(x = object)
    warning("Number of dimensions specified is greater than the number of cells. Setting dims to ", dims, " and continuing", immediate. = TRUE)
  }

  #(A7) 获取每个基因的 loading，就是组成PC时，每个基因的系数
  loadings <- Loadings(object = object[[reduction]], projected = FALSE)

  #(A8) 表示PC的基因名，一般是高变基因
  reduc.features <- rownames(x = loadings)
  # 如果基因数小于3个，则报错
  if (length(x = reduc.features) < 3) {
    stop("Too few features")
  }
  # 如果 基因数 x 每次抽样百分比 小于3，则警告
  if (length(x = reduc.features) * prop.freq < 3) {
    warning(
      "Number of variable genes given ",
      prop.freq,
      " as the prop.freq is low. Consider including more variable genes and/or increasing prop.freq. ",
      "Continuing with 3 genes in every random sampling."
    )
  }

  # (A9) 获取这些 表示PC的基因们对应的 scale.data
  data.use <- GetAssayData(object = object, assay = assay, slot = "scale.data")[reduc.features, ]





  #(A10) 读取日志: RunPCA.RNA
  
  # Seurat 的日之类一直没看 ...  能记录日期、命令本身、参数等信息。 //todo

  #> DefaultAssay(pbmc)
  #[1] "RNA"
  #> pbmc[["RunPCA.RNA"]]
  #Command: RunPCA(sce, features = VariableFeatures(object = sce))
  #Time: 2022-01-09 16:25:23
  #assay : RNA 
  #features : PPBP LYZ S100A9 IGLL5 ...
  #

  #> slotNames(pbmc@commands$RunPCA.RNA) #只有5个 slot，但是最后一个是list，可以记录很多参数
  #[1] "name"        "time.stamp"  "assay.used"  "call.string" "params"
  # > str(pbmc@commands$RunPCA.RNA@params)
  #List of 11
  #$ assay          : chr "RNA"
  #$ features       : chr [1:2000] "PPBP" "LYZ" "S100A9" "IGLL5" ...
  #$ npcs           : num 50
  #$ rev.pca        : logi FALSE
  #$ weight.by.var  : logi TRUE
  #$ verbose        : logi TRUE
  #$ ndims.print    : int [1:5] 1 2 3 4 5
  #$ nfeatures.print: num 30
  #$ reduction.name : chr "pca"
  #$ reduction.key  : chr "PC_"
  #$ seed.use       : num 42

  # 这个干啥的?
  rev.pca <- object[[paste0('RunPCA.', assay)]]$rev.pca
  #> pbmc[['RunPCA.RNA']]$rev.pca
  #[1] FALSE

  #是否乘以变异?
  weight.by.var <- object[[paste0('RunPCA.', assay)]]$weight.by.var
  # > pbmc[['RunPCA.RNA']]$weight.by.var
  #[1] TRUE





  #(A11) 开始对细胞抽样，做PCA，获取 基因的 系数矩阵
  fake.vals.raw <- my.lapply( #可以认为一个普通的 lapply
    X = 1:num.replicate, #重复抽样次数
    FUN = JackRandom, #抽样函数，返回的 loading 矩阵, 指定基因，指定PC[r1.use:r2.use]

    scaled.data = data.use, #以下是FUN的具体参数。
    prop.use = prop.freq, #抽样百分比
    r1.use = 1,
    r2.use = dims,
    rev.pca = rev.pca,
    weight.by.var = weight.by.var,
    maxit = maxit
  )


  #(A12) 整理抽样计算的结果
  # 这里面的返回值是什么？嵌套太多，没搞明白。

  fake.vals <- sapply(
    X = 1:dims, #遍历每一个PC
    FUN = function(x) { #x=PC编号
      # 又嵌入一层 lapply()
      return(as.numeric(x = unlist(x = lapply( #解开 lapply的list 后，变为数字，返回的是什么呢？
        X = 1:num.replicate, #遍历每一次抽样重复
        FUN = function(y) { #y=抽样次数的编号
          return(fake.vals.raw[[y]][, x]) #返回第y次抽样时，PC_x 那一列的值
        }
      ))))
    }
  )
  # 强转为矩阵
  fake.vals <- as.matrix(x = fake.vals)


  #(A13) 计算经验P值，返回的还是一个矩阵
  jackStraw.empP <- as.matrix(
    my.sapply(
      X = 1:dims, #遍历每一个PC
      FUN = function(x) { #x为PC编号
        return(unlist(x = lapply( #
          X = abs(loadings[, x]), #见上文(A7), 该PC的真实PC系数矩阵，取第x列，取绝对值
          # 对该列的每个基因的系数做遍历
          FUN = EmpiricalP, #这个函数太粗糙了吧？ 就是直接看 绝对值比x值大的占总数据个数的比例

          nullval = abs(fake.vals[,x]) #取第x列，取绝对值
        )))
      }
    )
  )
  # 加列名
  colnames(x = jackStraw.empP) <- paste0("PC", 1:ncol(x = jackStraw.empP))

  #(A14) 创建 JackStrawData 对象
  # pbmc@reductions$pca@jackstraw
  jackstraw.obj <- new(
    Class = "JackStrawData",
    empirical.p.values  = jackStraw.empP,
    fake.reduction.scores = fake.vals,
    empirical.p.values.full = matrix()
  )

  #(A15) 添加对象到 'pca' 中，pbmc@reductions$pca@jackstraw
  JS(object = object[[reduction]]) <- jackstraw.obj

  #(A16) 记录日志
  object <- LogSeuratCommand(object = object)

  return(object)
}







==> 多apply嵌套的返回值，测试如下:

fake.vals.raw=lapply(
  X=1:3,
  FUN=function(x){
    return( iris[1:x, 1:4])
  }
)
fake.vals.raw
输出:
[[1]]
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2

[[2]]
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2

[[3]]
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2




fake.vals <- sapply(
  X = 1:4, #遍历每一个PC
  FUN = function(x) { #x=PC编号
    # 又嵌入一层 lapply()
    return(as.numeric(x = unlist(x = lapply( #解开 lapply的list 后，变为数字，返回的是什么呢？
      X = 1:3, #遍历每一次抽样重复
      FUN = function(y) { #y=抽样次数的编号
        return(fake.vals.raw[[y]][, x]) #返回第y次抽样时，PC_x 那一列的值
      }
    ))))
  }
)
fake.vals <- as.matrix(x = fake.vals)
输出:
> fake.vals
     [,1] [,2] [,3] [,4]
[1,]  5.1  3.5  1.4  0.2
[2,]  5.1  3.5  1.4  0.2
[3,]  4.9  3.0  1.4  0.2
[4,]  5.1  3.5  1.4  0.2
[5,]  4.9  3.0  1.4  0.2
[6,]  4.7  3.2  1.3  0.2



loadings=iris[100:102,1:4]
loadings
输出:
> loadings
    Sepal.Length Sepal.Width Petal.Length Petal.Width
100          5.7         2.8          4.1         1.3
101          6.3         3.3          6.0         2.5
102          5.8         2.7          5.1         1.9



jackStraw.empP <- as.matrix(
  sapply(
    X = 1:4, #遍历每一个PC
    FUN = function(x) { #x为PC编号
      return(unlist(x = lapply( #
        X = abs(loadings[, x]), #见上文(A7), 该PC的真实PC系数矩阵，取第x列，取绝对值
        # 对该列的每个基因的系数做遍历
        FUN = Seurat:::EmpiricalP, #这个函数太粗糙了吧？ 就是直接看 绝对值比x值大的占总数据个数的比例

        nullval = abs(fake.vals[,x]) #取第x列，取绝对值
      )))
    }
  )
)
jackStraw.empP
输出:
> jackStraw.empP
     [,1] [,2] [,3] [,4]
[1,]    0  1.0    0    0
[2,]    0  0.5    0    0
[3,]    0  1.0    0    0
















() JS() 读写 jackStraw 数据
./seurat-object-4.0.4/R/dimreduc.R:257:"JS<-.DimReduc" <- function(object, slot = NULL, ..., value) {

参数类型:
> class(pbmc[["pca"]])
[1] "DimReduc"
attr(,"package")
[1] "SeuratObject"



#' @rdname JS
#' @export
#' @method JS<- DimReduc
#'
"JS<-.DimReduc" <- function(object, slot = NULL, ..., value) {
  CheckDots(...)
  if (inherits(x = value, what = 'JackStrawData')) {
  	# 干活的就是这一句，保存在 pbmc@reductions$pca@jackstraw
    slot(object = object, name = 'jackstraw') <- value
  } else if (is.null(x = NULL)) {
    stop("A slot must be specified")
  } else {
    JS(object = JS(object = object), slot = slot) <- value
  }
  # 然后返回对象
  return(object)
}



#' @rdname JS
#' @export
#' @method JS JackStrawData
#'
JS.JackStrawData <- function(object, slot, ...) {
  CheckDots(...)
  slot <- switch(
    EXPR = slot,
    'empirical' = 'empirical.p.values',
    'fake' = 'fake.reduction.scores',
    'full' = 'empirical.p.values.full',
    'overall' = 'overall.p.values',
    slot
  )
  return(slot(object = object, name = slot))
}












()EmpiricalP()
./seurat-4.1.0/R/dimensional_reduction.R:1944:EmpiricalP <- function(x, nullval) {

#internal
EmpiricalP <- function(x, nullval) {
  return(sum(nullval > x) / length(x = nullval))
}

# 经验P值 = 空值中 大于 x 的个数，除以空值长度 










() length()
./seurat-object-4.0.4/R/dimreduc.R:526:length.DimReduc <- function(x) {
#' @describeIn DimReduc-methods The number of dimensions for a \code{DimReduc}
#' object
#'
#' @return \code{length}: The number of dimensions
#'
#' @export
#' @method length DimReduc
#'
length.DimReduc <- function(x) {
  return(ncol(x = Embeddings(object = x)))
}


测试:
> Embeddings(object = pbmc_small)[1:2,1:3]
                      PC_1       PC_2       PC_3
ATGCCAGAACGACT -0.77403708 -0.8996461 -0.2493078
CATGGCCTGTGCAT -0.02602702 -0.3466795  0.6651668


对应的函数是:
> Loadings(object = pbmc_small)[1:2,1:3]
             PC_1       PC_2        PC_3
PPBP   0.33832535 0.04095778  0.02926261
IGLL5 -0.03504289 0.05815335 -0.29906272










() JackRandom() 主力函数
对基因抽样，打乱cell标签，做PCA，返回基因的loading矩阵
./seurat-4.1.0/R/dimensional_reduction.R:2206:JackRandom <- function(

#internal
#
JackRandom <- function(
  scaled.data,
  prop.use = 0.01, #默认抽样 1%
  r1.use = 1,
  r2.use = 5,
  seed.use = 1,
  rev.pca = FALSE, #默认值 FALSE，cell x gene matrix 的列，也就是 gene 做PCA
  weight.by.var = weight.by.var, #默认T, 对 cell embeddings 乘以 var 权重 
  # //todo 这里是bug: 参数默认值给一个未知变量?!
  maxit = 1000
) {
  # (A1)如果有种子，则设置
  if (!is.null(x = seed.use)) {
    set.seed(seed = seed.use)
  }

  # (A2) 打乱基因顺序，抽样 1%，也就是2000个高变基因 抽取 20个
  rand.genes <- sample(
    x = rownames(x = scaled.data),
    size = nrow(x = scaled.data) * prop.use
  )
  # 如果抽取基因数小于3个，则重新抽取，抽取3个
  # //todo 为什么不直接在 上一个函数的 size=中直接使用 ifelse()呢？
  # make sure that rand.genes is at least 3
  if (length(x = rand.genes) < 3) {
    rand.genes <- sample(x = rownames(x = scaled.data), size = 3)
  }

  #(A3) 传入的 scaled.data 赋值给新变量
  data.mod <- scaled.data

  #(A4) 新矩阵洗牌: 每一行独立进行打乱。一行是一个基因，就是每一行把 cell id 打乱了。
  data.mod[rand.genes, ] <- MatrixRowShuffle(x = scaled.data[rand.genes, ])

  #(A5) 运行PCA，
  temp.object <- RunPCA(
    object = data.mod, # scaled.data 抽样后，按行打乱
    assay = "temp", #指定 临时 assay
    npcs = r2.use, #使用的PC数，默认5个
    features = rownames(x = data.mod), #基因名列表
    rev.pca = rev.pca,  # 默认按列 - gene 做 PCA
    weight.by.var = weight.by.var, # 是否对cell embedding 乘以 var权重
    verbose = FALSE, #不显示进度条等
    maxit = maxit #求svd时最大迭代次数
  )

  # 返回基因的系数矩阵，指定基因，指定PC
  return(Loadings(temp.object)[rand.genes, r1.use:r2.use])
}








()MatrixRowShuffle()
./seurat-4.1.0/R/utilities.R:2075:MatrixRowShuffle <- function(x) {

见本文 3.2-3.3








() ScoreJackStraw()

pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
影响范围: pbmc@reductions$pca@jackstraw$overall.p.values
> head(pbmc@reductions$pca@jackstraw$overall.p.values)
     PC         Score
[1,]  1 9.123154e-161
[2,]  2 1.468402e-114

> tail(pbmc@reductions$pca@jackstraw$overall.p.values)
      PC      Score
...
[19,] 19 1.00000000
[20,] 20 0.24803564


#./seurat-4.1.0/R/generics.R:602:ScoreJackStraw <- function(object, ...) {

#' Compute Jackstraw scores significance.
#'
#' Significant PCs should show a p-value distribution that is
#' strongly skewed to the left compared to the null distribution.
#' The p-value for each PC is based on a proportion test comparing the number
#' of features with a p-value below a particular threshold (score.thresh), compared with the
#' proportion of features expected under a uniform distribution of p-values.
#'
#' @param object An object
#' @param ... Arguments passed to other methods
#'
#' @return Returns a Seurat object
#'
#' @author Omri Wurtzel
#' @seealso \code{\link{JackStrawPlot}}
#'
#' @rdname ScoreJackStraw
#' @export ScoreJackStraw
#'
ScoreJackStraw <- function(object, ...) {
  UseMethod(generic = 'ScoreJackStraw', object = object)
}




#./seurat-4.1.0/R/dimensional_reduction.R:1857:ScoreJackStraw.Seurat <- function(


#' @param reduction Reduction associated with JackStraw to score
#' @param do.plot Show plot. To return ggplot object, use \code{JackStrawPlot} after
#' running ScoreJackStraw.
#'
#' @seealso \code{\link{JackStrawPlot}}
#'
#' @rdname ScoreJackStraw
#' @concept dimensional_reduction
#' @export
#' @method ScoreJackStraw Seurat
#'
ScoreJackStraw.Seurat <- function(
  object,
  reduction = "pca",
  dims = 1:5,
  score.thresh = 1e-5,
  do.plot = FALSE,
  ...
) {
  
  # 抛给了 pca 类
  object[[reduction]] <- ScoreJackStraw(
    object = object[[reduction]],
    dims = dims,
    score.thresh = score.thresh,
    ...
  )

  if (do.plot) {
  	# 这个参数检查函数，以后看，貌似很复杂 //todo
    CheckDots(..., fxns = 'JackStrawPlot')
    # 根据代码低耦合的原则，这里最好不写这个画图功能。
    # 默认不走这个if，后面单独画图语句。
    suppressWarnings(expr = print(JackStrawPlot(
      object = object,
      reduction = reduction,
      dims = dims,
      ...
    )))
  }

  #记录日志，返回
  object <- LogSeuratCommand(object = object)
  return(object)
}





#./seurat-4.1.0/R/dimensional_reduction.R:1836:ScoreJackStraw.DimReduc <- function(object, dims = 1:5, score.thresh = 1e-5, ...) {

#' @rdname ScoreJackStraw
#' @concept dimensional_reduction
#' @export
#' @method ScoreJackStraw DimReduc
#'
ScoreJackStraw.DimReduc <- function(object, dims = 1:5, score.thresh = 1e-5, ...) {
  JS(object = object) <- ScoreJackStraw( #又抛给了 JackStrawData 类
    object = JS(object = object), # pbmc@reductions$pca@jackstraw
    dims = dims,
    score.thresh = score.thresh,
    ...
  )
  return(object)
}





由每个PC每个基因的经验p值，计算PC打分的主力函数。
#./seurat-4.1.0/R/dimensional_reduction.R:1796:ScoreJackStraw.JackStrawData <- function(

#' @param dims Which dimensions to examine
#' @param score.thresh Threshold to use for the proportion test of PC
#' significance (see Details)
#'
#' @importFrom stats prop.test
#'
#' @rdname ScoreJackStraw
#' @concept dimensional_reduction
#' @export
#' @method ScoreJackStraw JackStrawData
#'
ScoreJackStraw.JackStrawData <- function(
  object,
  dims = 1:5,
  score.thresh = 1e-5,
  ...
) {
  CheckDots(...)
  # 获取该slot的数据，是一个矩阵
  # pbmc@reductions$pca@jackstraw@empirical.p.values[1:10,1:5]
  pAll <- JS(object = object, slot = "empirical.p.values")

  # 只选取某几列
  pAll <- pAll[, dims, drop = FALSE]

  # 转为数据框
  pAll <- as.data.frame(pAll)

  #添加一列 基因名
  pAll$Contig <- rownames(x = pAll)
  
  # 打分数据框 初始化
  score.df <- NULL



  ##########################
  # 核心检验步骤: prop.test
  ##########################

  # 遍历PC
  for (i in dims) { #i 是每个PC的下标
  	# 百分比检验?
    pc.score <- suppressWarnings(prop.test(
      x = c(
        length(x = which(x = pAll[, i] <= score.thresh)), #小于某阈值的个数
        floor(x = nrow(x = pAll) * score.thresh) #总个数 * 阈值，取更小的整数
      ),
      n = c(nrow(pAll), nrow(pAll)) #两个一样的值
    )$p.val) #我认为这就是卡方检验

    # 如果 小于阈值的个数为0，则记p=1
    if (length(x = which(x = pAll[, i] <= score.thresh)) == 0) {
      pc.score <- 1
    }

    # 如果 score.df 空，则新建df，包含2列
    if (is.null(x = score.df)) {
      score.df <- data.frame(PC = paste0("PC", i), Score = pc.score)
    # 如果 score.df 已经有内容，则新增行
    } else {
      score.df <- rbind(score.df, data.frame(PC = paste0("PC", i), Score = pc.score))
    }
  }

  # 新增一列，算覆盖赋值，去掉了PC前缀，只要PC编号
  score.df$PC <- dims

  # 转为矩阵
  score.df <- as.matrix(score.df)

  # 保存到 pbmc@reductions$pca@jackstraw$overall.p.values
  JS(object = object, slot = 'overall') <- score.df
  return(object)
}











()
# JackStrawPlot(pbmc, dims = 1:15)
#./seurat-4.1.0/R/visualization.R:3782:JackStrawPlot <- function(

对于每个PC，画QQ图: 每个基因的p值，对比均匀分布。同时，也确定每个PC的总体p值。
显著的PC应该明显偏离null 分布，分布在黑色虚线的左侧。
每个PC的总体p值，是基于 proportion test 的，比较 p值低于某个阈值(score.thresh)的基因数，以及 基于均匀分布，低于该阈值的期望基因数。

#' JackStraw Plot
#'
#' Plots the results of the JackStraw analysis for PCA significance. For each
#' PC, plots a QQ-plot comparing the distribution of p-values for all genes
#' across each PC, compared with a uniform distribution. Also determines a
#' p-value for the overall significance of each PC (see Details).
#'
#' Significant PCs should show a p-value distribution (black curve) that is
#' strongly skewed to the left compared to the null distribution (dashed line)
#' The p-value for each PC is based on a proportion test comparing the number
#' of genes with a p-value below a particular threshold (score.thresh), compared with the
#' proportion of genes expected under a uniform distribution of p-values.
#'
#' @param object Seurat object
#' @param dims Dims to plot
#' @param cols Vector of colors, each color corresponds to an individual PC. This may also be a single character
#' or numeric value corresponding to a palette as specified by \code{\link[RColorBrewer]{brewer.pal.info}}.
#' By default, ggplot2 assigns colors. We also include a number of palettes from the pals package.
#' See \code{\link{DiscretePalette}} for details.
#' @param reduction reduction to pull jackstraw info from
#' @param xmax X-axis maximum on each QQ plot.
#' @param ymax Y-axis maximum on each QQ plot.
#'
#' @return A ggplot object
#'
#' @author Omri Wurtzel
#' @seealso \code{\link{ScoreJackStraw}}
#'
#' @importFrom stats qunif
#' @importFrom scales hue_pal
#' @importFrom ggplot2 ggplot aes_string stat_qq labs xlim ylim
#' coord_flip geom_abline guides guide_legend
#' @importFrom cowplot theme_cowplot
#'
#' @export
#' @concept visualization
#'
#' @examples
#' data("pbmc_small")
#' JackStrawPlot(object = pbmc_small)
#'
JackStrawPlot <- function(
  object,
  dims = 1:5,
  cols = NULL,
  reduction = 'pca',
  xmax = 0.1,
  ymax = 0.3
) {
  #(A1)获取的是 矩阵 obj@empirical.p.values
  pAll <- JS(object = object[[reduction]], slot = 'empirical')
  #如果最大维数超过 经验p值的列，报错
  if (max(dims) > ncol(x = pAll)) {
    stop("Max dimension is ", ncol(x = pAll))
  }
  # 按参数，传入 PC 列
  pAll <- pAll[, dims, drop = FALSE]
  # 转为数据框
  pAll <- as.data.frame(x = pAll)
  # 宽变长，使用的内部函数 Melt，自己写的才是最稳定的!
  # 吐槽一下: dplyr/tidyr 等动不动就改接口，频率太高，写包的都不敢用。
  # 比如这里，Seurat 包的作者只敢在使用示例中使用tidyverse包，而不敢在源码中使用。
  data.plot <- Melt(x = pAll) #Melt函数见 源码解析 19-3.11
  # 加列名
  colnames(x = data.plot) <- c("Contig", "PC", "Value")



  #(A2) 取 obj@overall.p.values
  score.df <- JS(object = object[[reduction]], slot = 'overall')

  #行数(PC数) 小于 要显示的PC数，报错
  if (nrow(x = score.df) < max(dims)) {
    stop("Jackstraw procedure not scored for all the provided dims. Please run ScoreJackStraw.")
  }

  # 按列，取子集 
  score.df <- score.df[dims, , drop = FALSE]

  # 如果行数为0，则报错
  if (nrow(x = score.df) == 0) {
    stop(paste0("JackStraw hasn't been scored. Please run ScoreJackStraw before plotting."))
  }


  #(A3) 为经验p值添加新列: 总体p值
  data.plot$PC.Score <- rep(
    x = paste0("PC ", score.df[ ,"PC"], ": ", sprintf("%1.3g", score.df[ ,"Score"])),
    each = length(x = unique(x = data.plot$Contig)) #每个x的个数等于独特基因的个数
  )
  #强转为因子
  data.plot$PC.Score <- factor(
    x = data.plot$PC.Score,
    levels = paste0("PC ", score.df[, "PC"], ": ", sprintf("%1.3g", score.df[, "Score"]))
  )


  #(A4) 如果颜色为空，则使用默认颜色。 
  if (is.null(x = cols)) {
    cols <- hue_pal()(length(x = dims))
  }
  # 如果颜色长度小于 PC 长度，报错
  if (length(x = cols) < length(x = dims)) {
    stop("Not enough colors for the number of dims selected")
  }

  #(A5) 画图
  # 输入参数 sample ，然后使用 stat_qq(distribution = qunif) 
  gp <- ggplot(data = data.plot, mapping = aes_string(sample = 'Value', color = 'PC.Score')) +
    stat_qq(distribution = qunif) +

    # 其他绘图参数
    labs(x = "Theoretical [runif(1000)]", y = "Empirical") +
    scale_color_manual(values = cols) +
    xlim(0, ymax) +
    ylim(0, xmax) +
    coord_flip() + #交换x、y坐标及标签
    geom_abline(intercept = 0, slope = 1, linetype = "dashed", na.rm = TRUE) + #添加斜线，斜率1，过原点
    guides(color = guide_legend(title = "PC: p-value")) + #修改图例标题
    theme_cowplot() #设置主题
  #(A6) 返回图
  return(gp)
}

